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UNSTEADY  SEPARATION  CHARACTERISTICS  OF  AIRFOILS 
OPERATING  UNDER  DYNAMIC  STALL  CONDITIONS 


by 

Wolfgang  Geifller,*  Lawrence  W.  Carr,**  Tuncer  Cebeci*** 


Abstract 


Within  the  scope  of  a  joint  effort  between  DFVLR,  NASA- Ames  and 
Cal.  State  University  Long  Beach,  represented  by  the  authors, 
unsteady  viscous/inviscid  interaction  phenomena^ hava^baan> inves¬ 
tigated  on  airfoils  operating  under  dynamic  stall  conditions.  It 
is  well  known  from  experiments  that,  in  the  upstroke  region  of  a 
sinusoidally  oscillating  airfoil,  the  flow  remains  attached  up 
to  incidences  considerably  larger  than  the  static  stall  angle. 
Reversed  flow  areas  develop  close  to  the  airfoil  surface  without 
boundary  layer  separation.  These  complicated  flow  phenomena  are 
investigated  in  the  present  study  on  the  basis  of  coupling  pro¬ 
cedures  between  a  time-dependent  inviscid  panel  method  and  2-a 
unsteady  boundary  layer  codes.  Two  strategies  are  pursued: 

1.  Coupling  of  inviscid  panel  method  with  boundary  layer 
code^—  direct  mode.  v.  4 

2.  Strong  coupling  of  inviscid  panel  method  with  boundary 
layer  code^—J. inverse  mode. 

The  main  features  of  the  unsteady  time-marching  panel  method  and 
boundary  layer  codes  are  discussed.  Emphasis  is  placed  on  the 
investigation  of  numerical  stability  and  the  phenomenon  of 
unsteady  separation.  Future  steps  of  the  ongoing  cooperation  are 
outlined.  o 
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1.  Introduction 


The  dynamic  stall  problem  is  closely  related  to  helicopter  aero¬ 
dynamics.  On  the  retreating  blade  of  a  helicopter  in  forward 
flight,  the  flow  may  separate  and  reattach,  leading  to  very  com¬ 
plicated  time-dependent  flow  phenomena.  In  the  present  study  this 
problem  is  simplified  in  the  usual  way  by  representing  the 
rotating  blade  by  a  nonrotating  characteristic  airfoil  section 
undergoing  pitching  oscillations  about  the  quarter-chord  axis. 
Extensive  windtunnel  experiments  [1]  have  shown  that  the  follow¬ 
ing  list  of  parameters  is  of  the  highest  importance  for  dynamic 
stall  characteristics: 
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airfoil  shape  Mach  number 

Reynolds  number  reduced  frequency 

mean  anqle,  amplitude  type  of  motion 

3-d  effects  tunnel  effects 

A  calculation  procedure  should  have  the  capability  to  take  most 
of  these  parameters  into  account.  This  has  been  verified  in  the 
present  methods  with  the  exceptions  ox  3-d  and  tunnel  effects. 

Flg.1  shows  measured  lift  and 
pitching  moment  distributions 

as  functions  of  incidence  THE  EVENTS  OF 

together  with  the  related  DYNAMIC  STALL  ON 

flow  events  (2].  The  lift  THE  NACA  0012  AIRFOIL 
curve  shows  a  considerable 
extension  beyond  the  static  , 

stall  angle  during  upstroke.  j  JT 

In  this  regime  the  boundary  w|  | 

layer  on  the  upper  surface  g 

still  remains  attached.  Hot  ^ 

film  measurements  and  flow  I 

visualizations  have  shown  z 

however  that,  at  these  inci¬ 
dences,  the  flow  inside  the 
boundary  layer  is  partially 

reversed.  This  reversed  flow 
area  is  located  adjacent  to  J 

the  wall  and  may  spread  over  t 

most  of  the  upper  surface.  | 

Lift  and  moment  curves  show  ° 

the  characteristic  hysteresis 
loops. 

In  the  present  study  the  up¬ 
stroke  region  will  be  invest¬ 
igated  on  the  basis  of  coup¬ 
ling  procedures  between  an 
unsteady  inviscid  panel  meth-  Fig.1 

od  and  time-dependent  bounda¬ 
ry  layer  codes.  Recent  investigations  with  unsteady  finite  dif¬ 
ference  boundary  layer  procedures  [3,4,5]  have  shown  that  calcu¬ 
lations  within  reversed  flow  areas  are  possible.  In  the  following 
discussion  the  limits  of  these  methods  due  to  numerical  stability 
(CFL  condition)  as  well  as  due  to  unsteady  separation  (MRS  cri¬ 
terion)  will  be  outlined.  Two  different  strategies  of  coupling 
procedures  between  inviscid  and  viscous  methods  will  be  applied: 

•  weak  coupling,  boundary  layer  calculation  in  the  direct  mode 

•  strong  coupling,  boundary  layer  calculation  in  the  inverse 
mode. 


Some  typical  results  of  these  calculations  will  be  given.  Further 
steps  are  necessary  and  planned  in  the  future  to  improve  tbe 
various  codes. 


2.  Unitea 


Fig. 2a  shows  ths  noman- 
elsturs  of  profils  and 
wake  for  ths  unstsady 
psnsl  msthod.  Both  pro- 
filss  and  waks  surfacss 
ars  subdivided  into 
panels,  each  of  which 
is  represented  by  a 
constant  source  and 
vorticity  distribution. 
Ths  vorticity  distribu¬ 
tion  is  extended  into 
ths  unsteady  wake.  The 
wake  geometry  is  sim¬ 
plified  by  a  straight 
line. 


Asm 


The  incidence  variation 


s  :  Surtax  coordhate  (Comm  point} 
s‘  :  Strfoce  coordinate  (Integration  w viable) 
r  :  Unit  tangential  vector  (in  s- direction) 

X  •'  Uni  normal  vector  (cuter  normal  drection) 


(1)  a(T)  *  a0*a1sin«#T 
with  T  *  t  ^  ;  u*  ■ 


Fig. 2a  Nomenclature 


is  subdivided  into  time  steps.  For  each  time  step  the  kinematic 
flow  condition  (zero  normal  velocity  on  the  profile  surface)  is 
fulfilled  and  the  strengths  of  the  source  and  vorticity  distrib¬ 
utions  ars  determined  by  solving  a  linear  system  of  equations. 
Ths  resulting  relative  velocity  vector  used  in  the  boundary  layer 
codes  as  an  outer  boundary  condition  is  obtained  by 


♦  v  ♦  v 
q  T 


Ue(x,T) 


with  vkin  88  the  prescribed 
kinematic  velocity  and  Vq 


and  vT  as  the  induced  ve¬ 


locities  of  the  source  and 
vorticity  system. 


The  unsteady  Bernoulli 
equation  determines  the 
pressure  coefficient 


(3)  Cp(TI  - 


with  <0(T)  as  the  velocity 
potential . 


By  means  of  the  unsteady 
Kutta  condition 


<4>si|;i=pu'T»  -  cPt'T'J 

(si  surface  coordinate) 


AS*  :  Length  o(  Httr  Qeamt 
AS,  :  length  at  Profie  ftnet 
Aj,  :  Change  at  Profile  tbrtdfy 
doing  Tern  step  AT 

Fig.  2b  Relation  between  profile 
and  wake  vorticity 


the  vorticity  loading  of  the  wake  element  leaving  the  trailing 
edge  at  time  T  is  specified.  Fig. 2b  describes  the  relationship 
between  profile  and  wake  vorticity.  After  each  time  step  a  new 
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wake  element  is  shed  at  the 
trailing  edge  which  then  moves 
downstream  with  U„. 

A  more  detailed  description  of 
the  method  is  given  in  [6].  As  a 
typical  result  of  the  unsteady 
panel  method,  lift  and  moment 
distributions  for  three  differ¬ 
ent  reduced  frequencies  are  pre¬ 
sented  in  Figs.3a*3c  for  the  NACA 
0012  airfoil.  Due  to  unsteady 
wake  effects,  hysteresis  loops 
develop  even  in  the  present 
inviscid  case:  for  small  values 
of  w*  (Fig. 3a)  the  lift  loop  is 
traversed  in  counterclockwise 
direction.  For  higher  values  of 
a*  (Fig.  3c)  this  trend  is 
reversed.  For  an  intermediate 
value  of  u*  (Fig. 3b),  a  hystere¬ 
sis  loop  for  the  lift  has  disap¬ 
peared  entirely.  The  influence 
of  w*  on  the  pitching  moment  is 
minor. 


Unsteady  inviscid  lift  and  moment  distributions 


3.  Unsteady  (Laminar/Turbulent)  Boundary  Layor:  Oiract  Mods  (GoiBlor) 


3.1  Boundary  layer  Calculation  Procaduro 


In  addition  to  tha  axtarnal  in- 
viscid  valocity  distribution 
Ue(X,T) ,  corresponding  initial 
conditions  are  necessary  to 
start  the  boundary  layer  calcu¬ 
lation.  Fig. 4  shows  as  a  result 
of  the  inviscid  method  that  the 
front  stagnation  point  is  mov¬ 
ing.  To  start  the  boundary  lay¬ 
er  calculation  at  the  stagna¬ 
tion  point  for  all  time  steps, 
the  unsteady  boundary  layer 
equations  are  expressed  in  a 
stagnation-point- fixed  frame  of 
reference  [5]  with  the  coordi¬ 
nate 


s*S*'rriM»7 

S-CL2 
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AMES-A01 


(S)  X, 


!(8  * 


SSP(T)) 


(sSp  i  time -dependent  location 
**  of  the  stagnation  point). 


Fig. 4  Displacement  of 
stagnation  point 


with  ±  referring  to  the  upper/lower  surface,  respectively. 

The  kinematic  velocity  of  the  stagnation  point  is  then 
3x7 

•  +  v  as  —  1 


(6)  0 


Ikin 


In  tha  new  system  the  unsteady  boundary  layer  equations  yield 


(7)  Tx7  +  3^  * 


3ui  3ui  3ui  31 

(8)  -yjr  ♦  u,  *  Vl  ^  •  -j 

with  «  v1  VWe  ;  n1  ■  Y1  VSe 
and  the  boundary  conditions 


'l  Tx7  * 


[(’  *  £rl 


n.  ,  u,  ■  u.  ♦  u, 

'  '  1  'kin 

For  tha  turbulent  boundary  layer  calculations  tha  algebraic  eddy 
viscosity  formulation  (c^  inEq.(8))  by  Cebeci/Smith  [7]  has  beer, 
used.  In  tha  present  investigation  the  transition  point  is  fixe*c 
directly  behind  the  (moving)  stagnation  point.  For  the  boundary 
layer  calculations  in  the  direct  mode  a  finite  difference  grir. 
of  the  Crank-Nicolson  type  is  used.  Fig. 5  shows  as  a  typica . 
result  contour  lines  of  vorticity  ( I )  and  velocity  (u)  inside  the 
boundary  layer  along  the  profile  upper  surface  at  a  specific 
instant  of  nondimens ional  time 

(10)  ?  «  .  1.15 

corresponding  to  the  instantaneous  incidence  a  =  19.85°  t. 


Fig. 5  Contour  lines  of 
equel  vorticity  1 
(left)  end  equel 
velocity  u  (right) 


The  celculetion  ceme  to  e 
breekdown  et  the  right-hend 
me r gins  of  the  figures  et 
ebout  xi  =  0.71.  This 
breekdown  occurred  due  to 
the  violetion  of  the  numer- 
icel  stebility  condition 
described  in  the  following 
section. 


3.2  Stebility  Con  side  rations 

As  long  ss  the  flow  inside 
the  boundery  leyer  is  non- 
reversed  the  Crenk-Nicolson 
finite  difference  scheme 
epplied  for  the  present  2 -d 
unsteedy  boundery  leyer 
equetions  is  unconditionel- 
ly  steble.  It  is  shown  in 
Fig. 6  however  thet  stebili¬ 
ty  limits  exist  in  regions 
where  the  flow  is  reversed. 
If  the  cherecteristics  of 
the  points  1,  2,  3  where  the 
boundery  leyer  quentities 
ere  known  ere  bent  beck- 
werds  such  thet  no  informe- 
tion  cen  reech  the  midpoint 
of  the  difference  molecule. 


Fig. 6  Stebility  bounds  of  cen¬ 
tral  difference  scheme 


the  calculation  braaka  down  (violation  of  CFL  condition).  It  is 
obvious  that  tha  nosh  also  ratic  Ax/AT  is  ths  limiting  valuo.  To 
•xcaod  ths  stability  bounds  ths  tins  stsp  AT  has  to  bs  rsducsd 
insids  rsvsrssd  flow  arsas. 


Fig. 7  Maximum  reversed  flow  velocity: 

investigation  of  the  CFL  condition 


Fig. 7  shows  maximum  reversed 
flow  velocities  inside  the 
boundary  layer  for  two  dif¬ 
ferent  numbers  of  sub-time 
steps  (original  time  step: 
AT  *  0.01).  Ten  sub-time 
steps  are  not  sufficient. 
The  boundary  layer  calcula¬ 
tion  breaks  down  after  a 
few  oscillations.  With 
twenty  sub- time  steps,  how¬ 
ever,  a  smooth  increase  of 
(u)  is  achieved. 

•  "V 


3.3  Unsteady  Separation 

If  the  numerical  stability 
bound  has  been  removed  due 
to  a  corresponding  reduc¬ 
tion  of  the  time  step,  the 
calculation  may  run  into  a 
limit  which  must  be  inter¬ 
preted  as  unsteady  separa¬ 
tion.  Fig. 8  shows  the  well- 
known  steady  separation 
criterion  and  the  more  com¬ 
plicated  unsteady  MRS 
(Moore,  Rott  and  Sears) 
criterion.  In  the  unsteady 
case,  separation  is  obtained  at  the  position  of  zero  shear  stress 
off  the  wall  when  the  maximum  reversed  flow  velocity  at  this 
position  reaches  the  separation  velocity.  The  MRS  criterion  can¬ 
not  be  checked  in  e  simple  manner  during  the  ongoing  calculation 
due  to  the  unknown  separation  velocity. 


V  Ska ay  Boundary  laytr  (2d) 
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Fig. 8  Separation  criteria 
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Fig. 9  Reversed  velocity  vectors  end  integral  curves 
(direct  node) 


Fig. 9  is  an  attempt  to  explain  the  MRS  criterion  by  means  of  the 
numerical  results  of  the  maximum  reversed  flow  velocities  inside 
the  boundary  layer.  The  upper  plots  show  vectors  of  the  maximum 
reversed  flow  velocities  for  two  different  amplitudes  of  oscil¬ 
lation.  The  lower  plots  display  the  corresponding  integral 
curves.  These  integral  curves  merge  from  both  sides  into  one  line 
forming  an  envelope  which  is  identified  as  the  unsteady  sepa¬ 
ration  line.  The  tangents  of  these  lines  correspond  to  the 
unsteady  separation  velocity.  Similar  behavior  has  already  been 
observed  in  3-d  steady  boundary  layer  calculations  [8]  for  the 
integrated  wall  shear  stress  vectors. 


Fig.  10  displays  a  comparison  of  calculated  reversed  flow  areas 
with  experimental  data  (1]  for  the  deep  dynamic  stall  regime. 
Reattachment  is  approximated  by  the  lines  tw  =  0.  For  the  higher 
Re-numbers  the  correspondence  between  the  calculations  and  the 
experimental  data  is  quite  sufficient.  Deviations  occur  for  the 
low  Re-number  case  where  the  assumption  of  a  fixed  transition 
point  is  no  longer  valid. 


-V»y  ■ 


Fig.  10  Separation  and  reattachment  loops:  deep  dynamic 
stall  conditions 


Finally,  Fig.  11  shows  how 
the  results  of  the  bounda¬ 
ry  layer  calculation  may 
be  used  for  a  modified 
panel  method  to  take  the 
location  as  well  as  the 
speed  of  unsteady  separa¬ 
tion  into  account.  A  cor¬ 
responding  modification  of 
the  unsteady  Kutta  condi¬ 
tion  for  two  vortical 
wakes  has  been  discussed 
by  Sears  [9).  The  intro¬ 
duction  of  two  separating 
wakes  in  the  unsteady  pan¬ 
el  method  is  straightfor¬ 
ward. 


4.  Unsteady  (laminar/Turbulent)  Boundary  Layer:  Inverse  Mode  (Cabaci/Carr) 
4.1  Boundary  Layer  Calculation  Procedure 

Unsteady  boundary  layer  calculations  extending  into  the  reversed 
flow  area  for  some  model  flows  have  been  presented  in  [3]  with 
the  application  of  Keller's  box  method.  Ref.  [10]  describes  ar. 
inverse  procedure  for  steady  separating  flows  over  profiles.  A 
combination  of  these  two  methods  has  been  developed.  The  velocity 
distribution  obtained  from  the  unsteady  panel  method  presented 
in  Section  2  serves  again  as  the  outer  boundary  condition. 


Change  of  Airfoil  Circulation 


Fig.  11  Modification  of  unsteady 
panel  method 


To  overcome  the  problem  of  tho  oacillacing  stagnation  point  tho 
boundary  laysr  calculation  now  starts  from  the  instantansous 
positions  of  ths  stagnation  point  on  a  quasi-steady  basis  until 
a  proscribed  position  is  reached,  whsrs  ths  unsteady  inverse 
calculation  starts  do%mstreaa.  The  unsteady  boundary  layer 
equations,  Eqs.(7),(8),  rsaain  unchanged  except  that  the  stagna¬ 
tion-point-  fixed  coordinate  x1  is  replaced  by  the  profile-fixed 
coordinate  x.  The  boundary  conditions  (Eq.(9))  are  simply  changed 
by  replacing  U1kin  with  sero. 

The  calculation  is  now  allowed  to  work  in  the  inverse  mode.  The 
external  velocity 


<U)  Uev(x,T) 
with 

(12)  6Ue(x,T) 


Ue(x,T)  ♦  6Ue(x,T) 


d 


(Ue  6*) 


dx' 


is  modified  by  means  of  the  blowing  velocity 
<13>  VB  *  2K  (ue  * 


Eq.(12)  is  the  Hilbert  Integral  representation  [10],  which  takes 
the  upstream  (xa)  and  downstream  (x^)  displacement  effects  of  the 
boundary  layer  at  the  present  position  x  into  account.  The  cal¬ 
culation  is  done  by  several  sweeps  along  the  profile  surface 
where  the  downstream  displacement  values  are  taken  from  the  pre¬ 
vious  sweep.  This  sweep  process  is  performed  for  each  time  step. 


The  standard  box  method  is  applied  in  regions  where  the  flow  is 
directed  downstream.  The  zigzag  scheme  is  used  in  regions  of 
reversed  flow  inside  the  boundary  layer.  Fig. 12  displays  the 
stability  bounds  of  the  zigzag  scheme  inside  the  reversed  flow 
area  in  correspondence  with  the  central  difference  scheme  of 
Fig. 6.  The  influence  region  is  now  twice  as  large  as  that  of  the 
central  scheme.  The  zigzag  scheme  requires  therefore  less  time 
steps  to  fulfill  the  CFL  condition  but  it  is  not  unconditionally 
stable  either. 


Fig. 12  Stability  bounds  of  zigzag  difference  scheme 


4.2  Coupling  with  the  Inviscid  Method 


The  strong  global  coupling  betvscn  ths  inviscid  panel  method  and 
the  boundary  layar  coda  is  carried  out  itarativaly:  inviscid 
valocitias  Ue(x,T)  ara  first  calculated  for  a  complata  cycle  of 
oscillations.  This  velocity  distribution  is  than  used  as  the 
outer  boundary  condition  for  the  boundary  layar  coda.  Due  to  the 
inverse  calculation  procedure  the  external  velocity  is  modified 
(Eqs. (11). (12) ) .  The  bloving  velocities  vB(x,T),  Eq.(13),  are 
then  taken  into  account  in  the  panel  method  modifying  the  kine¬ 
matic  flow  condition.  A  new  external  velocity  distribution  is 
calculated  now  including  the  effect  of  boundary  layer  displace- 
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Fig.  13  Maximum  reversed  velocity  vectors  and 
integral  curves  (inverse  mode) 


results  from  the  boundary 
layer  calculation  in  the 
direct  mode.  In  the  in¬ 
verse  case,  the  limiting 
process  of  forming  enve¬ 
lopes  of  the  integral 
curves  is  avoided.  The 
external  velocities  are 
now  allowed  to  adjust  to 
the  development  of  the 
boundary  layer  displace¬ 
ment  effect. 

Fig.  14  presents  integral 
curves  of  the  velocity 
vectors  inside  the  bound¬ 
ary  layer  and  along  the 
profile  upper  surface  for 
three  different  incidenc- 
ces  during  the  oscillato¬ 
ry  cycle.  These  curves 
are  not  to  be  interpreted 
as  streamlines.  They  are 
integrated  from  the  in¬ 
stantaneous  velocity 
(u,v)  distributions  in¬ 
side  the  boundary  layer. 
In  Fig. 14  only  a  very 
narrow  region  close  to 


•ei2 


Fig.  14  Recirculation  areas  inside 
the  boundary  layer 


the  wall  is  displayed: 

At  about  a  =  15°  upstroke 
(top  diagram)  a  recircu¬ 
lation  area  develops  near 
the  trailing  edge  which 
rapidly  increases  in  size 
up  to  the  maximum  inci¬ 
dence  of  a  =  18*  (middle 
plot).  Beyond  the  maximum 
incidence  the  recircula¬ 
tion  area  still  has  a 
large  spatial  extension 
down  to  «  =  15°  (bottom 
plot)  due  to  hysteresis 
effects.  The  plots  show 
further  that  the  integral  curves  are  bent  along  the  recirculation 
zone  and  compressed  toward  the  trailing  edge.  The  recirculation 
zone  stays  completely  on  the  airfoil.  Experimental  investigations 
have  shown,  however,  that  the  recirculation  area  is  extended 
beyond  the  trailing  edge.  This  deficiency  of  the  numerical  method 
is  mainly  due  to  the  neglection  of  the  wake  displacement  which 
is  a  sink  effect  on  the  airfoil  influencing  the  Kutta  condition. 
Future  extension  of  the  present  method  should  therefore  suffi¬ 
ciently  include  the  wake  displacement. 

Finally,  Fig. IS  displays  the  calculated  lift  and  moment  curves  for 
a  moderate  incidence  case  compared  with  experimental  data  (11). 
The  dotted  slopes  indicate  the  results  of  the  inviscid  panel 
method.  The  solid  curves  show  the  results  including  the  inverse 
boundary  layer  calculation  after  one  global  cycle.  For  this  case 
where  no  reversed  flow  occurs  during  the  oscillatory  cycle,  the 


agreement  between  theory  end  date  ia  very  good.  It  is  assumed 
that  correspondingly  good  results  will  be  achieved  at  higher 
incidence  variations  if  unsteady  wake  effects  are  taken  into 
account  in  the  coupling  procedure. 


Fig. IS  Unsteady  lift/moment  distributions: 

inviscid,  viscous,  experimental  data 


5.  Concluding  Remarks 

For  the  investigation  of  unsteady  viscous/lnviscid  interaction 
phenomena  coupling  procedures  between  a  panel  method  and  unsteady 
boundary  layer  codes  have  been  developed.  Working  in  the  direct 
mode  where  the  outer  inviscid  velocity  distribution  is  unchanged, 
it  has  been  shown  that  the  calculation  can  be  extended  into 
regions  of  reversed  flow.  The  asymptotic  behavior  of  boundary 
layer  quantities  within  the  reversed  flow  area  signals  unsteady 
separation.  By  investigation  of  the  maximum  reversed  flow  veloc¬ 
ities  it  has  been  shown  that  the  corresponding  integral  curves 
form  envelopes  which  are  interpreted  as  the  locations  of  unsteady 
separation.  The  information  on  separation  locations  and  sepa¬ 
ration  velocities  may  serve  as  inputs  for  a  modified  panel  meth¬ 
od. 

Working  in  the  inverse  mode,  the  external  velocity  distribution 
is  adjusted  to  the  development  of  the  boundary  layer  displacement 
thickness.  Analogous  integral  curves  of  the  maximum  reversed  flow 
velocities  do  not  show  the  asymptotic  behavior  of  the  direct 
mode.  A  singularity  in  the  boundary  layer  solution  is  obviously 
avoided  in  the  inverse  mode.  For  moderate  incidence  variations 
the  lift  and  moment  curves  agree  well  with  experimental  data.  For 
higher  incidences  the  effect  of  the  wake  has  a  severe  influence 
on  the  results.  Modification  of  the  method  to  consider  the  boun¬ 
dary  layer  along  the  unsteady  wake  as  well  as  a  representation 
of  the  wake  displacement  by  a  sink  distribution  in  the  panel 
method  will  be  necessary  future  activities. 
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